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1.  Introduction 


In  1880  Pierre  and  Jacques  Curie  experimentally  demonstrated  the  direct  piezoelectricity  effect  in 
which  a  mechanical  load  generates  an  electric  charge.  In  1881  Lippman  (7),  based  on 
thermodynamic  considerations,  postulated  the  inverse  effect,  where  an  electric  field  generates  a 
mechanical  response.  His  prediction  was  subsequently  confirmed  by  the  Curies.  Voigt  presented 
the  first  general  thermodynamic  theory  of  piezoelectricity  during  the  1890-1894  time  frame  (2). 
Since  then,  various  theories  of  linear  piezoelectricity  have  arisen  based  upon  assumptions 
concerning  both  deformation  and  the  material  constitutive  response.  The  majority  of  these 
theories  begin  with  linear  elasticity  and  the  assumption  of  a  symmetric  stress  tensor  (2).  Some 
researchers  have  tried  to  reformulate  the  governing  equations  in  terms  of  a  Cosserat  continuum 
(4).  Additional  simplifying  assumptions  are  subsequently  used  to  develop  equations  relevant  for 
plates,  shells,  and  beams  (5-7). 

The  behavior  of  piezoelectric  materials  in  non- structural  applications  has  been  investigated 
extensively.  However,  these  investigations  almost  exclusively  employ  linear  constitutive 
relations.  For  example,  piezoelectricity  is  undergoing  a  resurgence  in  both  fundamental  research 
and  technical  applications  (8-12).  However,  investigations  in  the  nonlinear  basic  theory  for 
piezoelectricity  have  been  limited.  Nelson  (13),  Toupin  (14),  and  Tiersten  (15,  16)  studied  the 
nonlinear  theory  of  dielectrics.  Norwood  et  al.  (17)  and  Kulkami  and  Hanagud  (18)  used  a 
Neo-Hookean  constitutive  relation  to  model  the  response  of  piezoelectric  ceramics.  Pai  et  al. 

(19)  considered  the  dependence  of  the  piezoelectric  strain  parameters  upon  the  strain  in 
formulating  a  plate  theory  of  piezoelectric  laminates.  Joshi  (20)  considered  the  nonlinear 
constitutive  relations  for  piezoelectric  materials,  where  a  concise  expression  was  given.  Tiersten 
(21)  investigated  the  nonlinear  problems  of  thin  plates  subjected  to  large  driving  voltages. 
Recently,  Patel  et  al.  used  finite  element  techniques  to  solve  a  particular  form  of  the  nonlinear 
piezoelectric  equations  with  a  nonlinear  stress-strain  relationship  but  linear  electrostatics  (22). 
However,  they  assume  that  the  effects  of  the  nonlinear  constitutive  can  be  neglected  and  only 
include  nonlinear  strain  effects.  Based  on  the  theory  of  invariants,  from  invariant  polynomial 
constitutive  relations,  Yang  and  Batra  (23)  investigated  the  second-order  theory  for  piezoelectric 
materials  with  symmetry  class  6mm  and  class  mm2.  Feng  et  al.  (24),  using  results  from  Kiral 
and  Eringen  (25),  developed  the  relations  for  symmetry  classes  6mm  and  3m  including  both 
nonlinear  stress-strain  and  electrostatic  constitutive  relations. 
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In  this  report,  we  present  the  Feng  et  al.  (24)  equations  in  a  form  suitable  for  numerical  solution. 
These  equations  are  then  specialized  for  one  dimensional  (1-D)  use.  Numerical  solutions  are 
verified  by  comparison  with  exact  solutions  of  the  linear  piezoelectric  equations  obtained  using 
Laplace  transform  techniques.  Finally,  the  response  of  the  full  nonlinear  equations  to  both  step 
pressure  and  step  voltage  boundary  conditions  are  examined. 


2.  Finite  Deformation  Piezoelectricity 


The  aforementioned  theories  are  typically  obtained  by  formulating  the  equations  of 
piezoelectricity  starting  with  principles  from  continuum  mechanics.  This  approach  provides  a 
general  formulation  with  a  clearer  understanding  of  the  restrictions  on  deformation  and  material 
constitutive  response  imposed  upon  the  resulting  governing  equations.  The  equations  of  motion, 
in  the  current  (Eulerian)  configuration,  can  be  written  as  (26) 

Tki,k  +  ph  =  pui  ,  (1) 

where  p  is  the  material  density,  bi  is  the  body  force,  ip  is  the  particle  acceleration,  and  T y  is  the 
total  stress.  Lowercase  Latin  subscripts  indicate  quantities  that  refer  to  the  current  configuration 
while  uppercase  Latin  subscripts  indicate  quantities  that  refer  to  the  reference  configuration  with 
subscripts  ranging  from  1-3.  Partial  differentiation  is  denoted  by  a  comma  before  the  index.  In 
equation  1,  symmetry  of  the  stress  tensor  is  not  required  although  we  include  it  in  the  following. 
Additionally,  there  are  no  assumptions  concerning  the  form  of  the  constitutive  relationships. 
Including  mechanical  and  electromagnetic  forces,  the  total  stress,  Tki,  can  be  written  as  (23) 

Tki  =  T'f ^  +  e0EkEl  —  -e0EmEmSki  ■  (2) 

Tc  is  the  Cauchy  stress  and  the  remaining  terms  are  the  Maxwell  stress  expressed  in  the  current 
configuration  (27),  where  eo  is  the  electric  permittivity  of  free  space,  E  is  the  electric  field,  and  5ki 
is  the  Kronecker  delta. 

The  equations  of  motion  in  the  reference  (Lagrangian)  configuration  are  obtained  from  equation 
1 .  The  equivalence  of  the  divergence  in  different  reference  frames  is  given  by 

( Gki),k  —  -j(J  E k\G  ki)  ,k  ,  (3) 

where  G  is  an  arbitrary  second  rank  tensor,  FkK  =  xk,K  is  the  deformation  gradient,  and 
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J  =  det(F).  Cauchy  stress,  T^,  is  written  in  terms  of  the  second  Piola-Kirchoff  stress,  TKL  as 


Tkt  —  —Xk^KL^i^L  —  ~jFkKTKLFiL 


The  electric  field  and  electric  displacement  transformations  are  given  by  (28) 


Ek  —  XKkEK  —  Fk\Ek 


(4) 


(5) 


F>k  =  ~jXk,K  F)  k  =  ~jFkKDK  (6) 

where  E  and  D  are  the  electric  field  and  the  electric  displacement  in  the  reference  configuration, 
respectively,  and  E  and  D  are  the  fields  in  the  current  configuration. 


The  appropriate  transformation  for  the  polarization,  Pk,  has  been  the  subject  of  debate  (29,  30). 
Several  commonly  used  transformation  are  shown  in  table  1.  Yang  and  Batra  (31)  used  a  set  of 
transformations  that  are  very  different  from  the  traditionally  accepted  forms,  including  the 
transformation  for  the  electric  field.  Dorfmann  and  Ogden  (29)  used  a  form  that  is  more 
commonly  accepted.  Their  rationale  for  the  form  of  the  transformation  for  the  polarization  is  that 
the  form  can  be  selected  so  as  to  simplify  mathematical  manipulation  of  the  equations.  Clayton 
(30)  selected  a  form  of  the  transformation  for  the  polarization  based  on  this  premise.  As  shown 
by  Lax  and  Nelson  (28),  the  forms  of  the  transformations  for  all  the  field  variables,  including 
polarization,  are  not  arbitrary.  The  appropriate  transformation  for  the  polarization  is  determined 
by  consideration  of  conservation  of  charge,  which  leads  to 


Pk 


1 

J 


Xk.K^K  ~ 


1 

J 


Ek  k  n  k 


(7) 


where  P  and  II  are  the  polarization  in  the  current  and  reference  configurations,  respectively. 


Table  1.  Transformation  relations. 


Author 

Electric 

Electric 

Polarization 

Displacement 

Field 

Yang/Batra  (31) 

<C 

II 

Cl 

£ 

II 

Ha  =  JF^Pa 

Dorfmann/Ogden  (29) 

<c 

L ^ 

E^ 

II 

c 

Ea  =  FaAEa 

Ha  =  JF^Pa 

Clayton  (30) 

<c 

E^ 

II 

c 

Ea  =  FaAEa 

IIa  =  FaAPa 

Lax/Nelson  (28) 

1  ^ 

Eti 

h-i 

II 

c 

Ea  =  FaAEa 

Ha  =  JF^Pa 
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Using  the  relations  FK\:  =  XKk,  Jp  =  p0,  where  p0  is  the  density  referred  to  the  reference 
configuration  with  equations  3-7,  the  equation  of  motion  in  the  Lagrangrian  frame  can  finally  be 
written 

{TKLxiiL)tK  +  (JFKl(e0F^kF^  -  +  Pok  =  PoUi  ■  (8) 

Equation  8  represents  the  balance  of  forces  acting  on  a  volume  of  material.  There  are  12 
unknowns,  6  stress  components,  3  electric  field  components,  and  3  displacements.  The 
additional  equations  required  for  a  consistent  formulation  are  obtained  from  the  governing 
equations  relating  the  electric  field  components  and  the  piezoelectric  constitutive  laws. 

The  governing  equations  relating  the  electric  field  components  are  obtained  by  starting  with 
Gauss’s  law  with  no  free  charges, 

(Dk),k  =  0  ,  (9) 

where  D  is  the  electric  displacement.  Then,  using  the  constitutive  relationship  Dk  =  e{)Ek  +  Pk, 
equation  9  becomes 

(to JfkIfm\em  +  n K)tK  =  o  .  (io) 

The  Euler-Piola-Jacobi  identity  (32),  {  JFKxk)  K  =  0,  is  used  to  simplify  equation  10  to 

(Rk),k  +  eoJFKl(FM\EM),K  =  0  .  (11) 

Equations  8  and  1 1  are  the  equations  of  motion  and  Gauss’s  law  in  the  reference  configuration. 
Closed-form  solutions  to  the  full  equations  given  by  equations  8  and  1 1  do  not  appear  feasible. 

A  standard  approach  for  solving  these  equations  is  the  use  of  numerical  methods.  For  instance, 
the  Galerkin  method  is  often  applied  in  finite  element  methods.  Another  approach  is  to  reduce 
the  equations  to  one  dimension  and  seek  exact  solutions.  As  shown  below,  we  will  illustrate  both 
of  the  approaches.  Accordingly,  we  first  express  the  governing  equations,  8  and  11,  in  their 
corresponding  weak  forms.  Then,  we  reduce  the  equations  to  their  1-D  representations.  Finally, 
we  can  also  obtain  solutions  to  these  equations  in  one  dimension  using  Laplace  transform 
techniques,  assuming  that  the  displacements  are  infinitesimal  and  that  the  piezoelectric 
constitutive  relations  are  linear.  We  then  can  compare  solutions  to  the  full  equations  expressed  in 
the  weak  form  with  the  1-D  exact  solutions. 


4 


3.  Weak  Form  Expressions 


The  weak  form  of  the  governing  equations  is  obtained  using  Galerkin’s  method  (33).  Applying 
this  technique  to  equations  8  and  1 1  leads  to  (34) 

—  [  (TklXi.l  +  {J F K\(eQF MkF N]  —  -e0FM1mFN^n5ki)EMEN))vitKdn.0  +  j  p0biVidU0 

J  Qo  ^  ^  Ho 

—  [  poUiVidfl0+  j  (TklXi.l  +  (JFj^eoF^F^  —  -e0FM1mFN^n5ki)EMEjy))riKVidr0 

Jn  o  J  r0  1 

—  [  (n^c  +  eoJFKl(FM1kEM))^,Kdflo  +  f  (11*  +  e0JFKl(FM\EM))nK^dY0  =  0  (12) 

Jfio  Jr  o 

where  Vi  and  <f>  are  arbitrary  displacement  and  scalar  electric  potential  test  functions,  respectively, 
f20  is  the  domain  of  the  reference  configuration  with  boundary  T0,  and  nK  is  the  outward  normal 
on  the  boundary.  Assuming  the  response  is  magnetostatic,  the  electric  field,  E,  can  be  expressed 
as  the  gradient  of  <I> 

E  =  -V<f>  .  (13) 

The  variational  statement  can  therefore  be  written 

—  j  {Tkl%i,l  +  {JFKl(eoFM\FNj  —  -e0FM1mFNln5ki)^,M^,N))vi,KdFl.o  +  f  pobit’idU0 

—  [  poUiVidtt0+  [  (TklXi:l  +  (JFj^^eoFj^F^  -e0FM1rnFNln8ki)^,M^,N))'nKVidr0 

J  n0  J  To  ^ 

—  [  (nx  —  eo JFKl(FM1k$tM))®,Kd£lo  +  f  (II*  —  eoJFKl(FM1k^)M))nK^dTQ  =  0  . 

J  Uo  Jr0 


For  small  electric  fields,  equation  14  reduces  to 


/  TKLxitLVij<dFl0  +  /  p0biVidfl0  —  /  p0uiVidU  0  —  /  HK&Kdtt0 
'  n0  j  n0  j  n0  j  n0 

+  /  Tklx^ jj n k Vi dT o  +  /  UxTiK^dEQ  =  0  .  (15) 

J  r0  Jr0 


The  solution  of  these  equations  is  obtained  by  minimizing  each  equation  with  respect  to  each  test 
function. 
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To  complete  the  equations,  TKL  and  UK  are  related  to  the  strains  and  the  electric  field  by 
appropriate  constitutive  relations.  These  relations  are  obtained  from 


Tkl 

n  k 


<9E 


dr  kl 

«9E 


dEK 


(16) 

(17) 


where  E  is  the  free  energy  density  function.  Expressions  for  E  are  lengthy.  Feng  et  al.  (35), 
and  Kiral  and  Eringen  (25)  provide  explicit  expressions.  In  one  dimension,  these  expressions 
reduce  to 


cET  -  eE  +  ^CT2  -  l-QE2  -  gTE  , 

(18) 

eT  +  eE  +  ^T2  +  ^ qE 2  +  QTE 

(19) 

where  cE  is  the  second-order  elastic  stiffness  coefficient,  e  is  the  piezoelectric  coupling 
coefficient,  e  is  the  dielectric  permittivity,  C  is  the  third-order  elastic  coefficient,  0  is 
electro strictive  coefficient,  g  is  the  third-order  piezoelectric  coupling  coefficient,  and  q  is  the 
third-order  dielectric  permittivity.  T  is  the  finite  strain  measure.  The  components  of  T  are  given 
by 


r /j  = 


1  .  du j  duj 

2  KdXj  dXj 


dux  9uk 
'dXI~dX~J 


(20) 


which  in  one  dimension  reduces  to 


du  1  . du  .  9 

- b  -( - ) 

dZ  2ydZJ 


(21) 


Consequently,  with  deformation  restricted  to  one  dimension  and  neglecting  the  Maxwell  stress 
and  body  force  terms,  the  final  weak  form  expression  is  given  by 

(-T(l  +u,z)v,z  -Tl$,z+pouv)dX  +Tx,z  v  |[,  +11$  \l0  =0  .  (22) 

Numerical  solutions  of  equation  22  were  obtained  using  the  COMSOL  Multiphysics  analysis 
software  (36).  The  generalized  a  method  was  employed  as  the  time  dependent  solver.  As  this 
solver  is  only  A-stable,  it  exhibits  spurious  high  frequency  ringing  when  subjected  to  a  step 
loading  (37).  To  address  this  issue,  artificial  Rayleigh  damping  was  incorporated  by  adding  an 
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additional  weak  term  in  the  form  (38) 


,  f  Fdv  d  du. 

d«  =  -J  VRC 

to  the  left-hand  side  of  equation  22,  where  vR  is  an  adjustable  parameter  used  to  minimize  the 
oscillations  in  the  solution. 


4.  Laplace  Transform  Solutions 


Solutions  of  equation  22  can  be  verified  by  comparison  with  solutions  to  the  strong  form  of  the 
equations  that  are  either  exact  solutions  or  solutions  obtained  by  a  different  numerical  technique. 
When  considering  the  deformations  in  equations  8  and  1 1  to  be  infinitesimal  and  the  electric 
fields  to  be  small,  and  neglecting  body  forces,  it  is  possible  to  obtain  solutions  that  can  be 
considered  exact  for  all  practical  purposes  using  Laplace  transform  techniques.  With  the 
assumptions  given,  equation  8  reduces  to  the  standard  wave  equation 


d2u  dax 
^  dt2  dx 


(24) 


In  equation  24  and  the  following,  we  express  all  equations  in  the  more  commonly  found  partial 
derivative  forms  rather  than  indicial  notation.  The  1-D  constitutive  equations  for  a  linear 
piezoelectric  material  can  be  written  as  (39) 


(T  x  c  hDx , 
dx 


du 


(25) 


E x  =  -h—  +  Dx/e,  (26) 

dx 

where  crx  is  the  stress  in  newtons/m2,  cD  is  Young’s  modulus  (at  constant  electric  displacement)  in 
newtons/m2,  u  is  the  particle  displacement  in  meters,  Dx  is  the  electric  flux  density  in  the 
x-direction  in  coul/m2,  e  is  the  permittivity  in  farads/m,  h  is  a  piezoelectric  constant  in  V/m,  and 
Ex  is  the  x-component  of  the  electric  field  in  V/m  (N/coul).  This  form  of  the  constitutive 
relationship  is  more  amenable  to  analytic  solutions  than  other  equivalent  forms  of  the 
piezoelectric  constitutive  relationships.  However,  the  use  of  different  constitutive  relations 
requires  determination  of  the  appropriate  transformations  relating  the  constitutive  coefficients.  In 
equations  25  and  26,  cD  =  cE  +  eh  and  h  =  e/e.  These  transformations  are  discussed  in 
appendix  B. 
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Substitution  of  equation  25  into  equation  24  gives, 

d2u  cD  d2u  dDx 
dt 2  p  dx 2  ^  dx 


(27) 


Assuming  plane  wave  propagation,  and  that  there  is  no  free  charge  inside  the  piezoelectric 
medium.  Gauss’s  Law,  given  by  equation  9,  reduces  to 


dDx 

dx 


0, 


(28) 


Consequently,  equation  27  again  reduces  to  the  1-D  wave  equation  for  isotropic  elastic 
(nonpiezoelectric)  media  and  can  be  written  as 

d2u  cD  d2u 
dt2  p  dx2' 


(29) 


4.1  1-D  Plane  Wave  Solutions 

In  this  section,  we  find  solutions  to  equation  29  with  linear  constitutive  equations  given  by 
equations  25  and  26  for  the  boundary  conditions  listed  in  table  2.  Figure  1  depicts  the  location 
and  orientation  of  the  boundary  conditions.  The  solutions  are  obtained  by  assuming  the 
displacements  u  in  the  piezoelectric  medium  can  be  expressed  as  D’Alembert  functions,  e.g., 

u(x,  t)  —  F2(t  —  x/c)  +  Fi(t  +  x/c)  (30) 


where  c  =  \JcD/p. 


Table  2.  Boundary  conditions. 


Location 

Step  Voltage 

Resonance 

Step  Pressure 

x  =  0 

4>  =  H(t),T  =  0 

<f>  =  0,  T  =  sin{u>t ) 

ct>  =  0,T  =  ToH(t) 

x  =  l 

II 

o 

4 

II 

o 

-Q- 

II 

o 

II 

o 

II 

o 

e 

II 

o 

Using  the  fact  that  C[f{t  —  r)H(t  —  r)]  =  e~STF(s )  t  >  0,  Laplace  transform  of  equation  30 

gives 

jC[u(x,  t)}  =  u(x,  s)  =  e~ F-2{s)  +  e~ Fi(s)  .  (31) 
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Figure  1 .  1  -D  boundary  value  problem. 


Substitution  of  the  spatial  derivative  of  equation  30  into  equation  25  yields, 

cD 

ax  +  hDx  = -[-F^t  -  x/c)  +  F[(t  +  x/c)\  . 

To  find  the  potential  across  a  plate  of  thickness,  l,  integrate  equation  26: 

i 

V{t)  =  J  Ex(x,t)dx, 

o 


(32) 


(33) 


or 

V(t)  =  - h[u(l )  -  u( 0)]  +  lDx(t) /e,  (34) 

where  the  quantity  u(l)  —  w(0)  represents  the  relative  motion  of  the  plate  surfaces.  Letting 
ox  =  a(x,  t ),  Dx(t)  =  D(t )  and  taking  the  Laplace  transform  of  equation  32  gives 

a(x,  s)  +  h,D  =  —  [— F2(s)  +  se Fi(s)].  (35) 

Laplace  transform  of  the  stress-free  boundary  at  0,  d(0,  s)  —  0  results  in 

cD 

hD  =  -~[-sF2{s)  +  sF1(s )]  (36) 

and  of  the  stress-free  boundary  at  /,  a (l,  s)  =  0  results  in 

hD  =  —  [— se~  F2(s)  +  se^Fi(s)].  (37) 
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Simultaneous  solution  of  equations  36  and  37  gives 


F2(s)  =  —e  <■  /•’,(*) 


(38) 


Taking  the  Laplace  transform  of  equation  34  and  substitution  of  equation  31  gives 

V(s)  =  Mi[e^F2(s)  +  e  c  Fi(s)  -  F2(s)  -  F^s)}  +  ID/e, 

or 

V'(.s)  =  /i[F2(s)(1  -  6#)  +  /'i(.s)(l  -  e« )]  +  ID/e. 

Substitution  of  equation  38  into  equation  36  yields 

-chDe^ 


(39) 


(40) 


or 


F,s)  = 


FiM  = 


cDs(  1  +  e  c ) 
chD 


cDs(  1  +  e  = ) 

Substitution  of  equations  41  and  42  into  equation  40  results  in  an  expression  for  the  Laplace 
transformed  potential: 

l  2ch2  tanh  (|^) 


(41) 


(42) 


y(s)  =  d 


Ms 


(43) 


For  a  step  voltage  input  applied  to  a  stress  free  plate,  V (t)  =  V0H(t )  and  V ( s )  =  V^/s.  Also 
J(f)  =  dQ(t)/dt  =  dD(t)/dt,  so  I  =  sQ  =  sD.  The  Laplace  transformed  current  can  be  written 
as 


I(s)  =  sQ(s)  =  sD(s )  = 


scDe  Vo 


cDsl  —  2ch2e  tanh 


(44) 


Also,  from  equation  35,  we  have  on  rearrangement 


cr(x,s)  =  —  [— se  c  F-2(s)  +  se  c  Fi(s)]  —  h.D. 


(45) 


Substitution  of  equations  41  through  43  into  equation  45  results  in  the  Laplace  transformed 
stresses: 

2hcDe sinh  (£^1)  sinh  (f )  C0 

<J(X,S)  =  - 7V7 - — - 77 — r-.  (46) 

V  '  2c/i2esinh  (^)  —  cDsl  cosh  (^) 

The  time  domain  solutions  are  found  by  numerical  inversion  of  the  Laplace  transform  using  the 
Dubner-Abate-Crump  (DAC)  algorithm  described  by  Crump  (40).  The  effects  of  Gibbs’ 
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phenomena  are  mitigated  using  the  so-called  Lanczos  a- factors  with  2048  terms  and  a  tolerance 
equal  to  10-4  (43).  Following  similar  procedures,  the  Laplace  transformed  stresses  for  the 
response  to  a  sinusoidal  voltage  and  a  step  pressure  are  given  by 


(2sujhcD  eV0sinh(  )sinh( (ffy)) 

^  ^  (s2  +  cu2)(— lcD  scosh(j^))  +  (2ch2esinh(j^)) 

and 

(lcD  scosh)^)  —  ch2esinh(l-^))Po 
lcDs2cosh( —  ch2sesinh(l-f) 

4.2  Plane  Waves  in  a  Free-free  Bi-material  Plate 


(47) 


(48) 


Finally,  we  can  also  obtain  solutions  for  a  bi-material  plate  (figure  2)  using  Laplace  transforms. 
Assuming  that  the  displacements  in  materials  1  and  2  are  given  by  u  \  and  u2,  respectively,  we  can 
write 


ui(x,t)  =  F2(t  -  x/ci)  +  Fi(t  +  x/ci ) 
u2(x,t)  =  F4(t  -  x/c2)  +  F3(t  +  x/c2). 


(49) 


O 


-  JC 


Figure  2.  Two-layer  boundary  value  problem. 
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Substitution  of  the  spatial  derivatives  of  equation  49  into  equation  25  yields, 


ai(x,t)  =  cf  ^  -  /ii-Di  =  -  x/ci]  +  F[(t  +  x/ci]  -  h1D1 

(r2(x,t)  =  cf  ^  -  h2D2  =  -  x/c2]  +  Fg(f  +  x/c2]  -  /i2£>2- 


(50) 


Taking  Laplace  transforms  of  the  stress  and  displacement  boundary  conditions  with  reference  to 
the  coordinate  system  shown  in  figure  2,  together  with  the  assumption  that  there  are  no  free 
charges  in  the  dielectrics,  i.e.,  (^Di  —  D 2j  •  h  —  0,  and  D\{s)  =  D2(s)  =  D(s),  permits 
determination  of  the  unknown  functions  F\ .  F-> ,  F:> ,  and  F4,  in  equations  49  via 


<xi(l  i,s)  =  <r2(li,s) 


ui(h,s)  =  u2(h,s) 


cri(0,  s)  =  0 


(51) 


or 


cr2(l,s)  =  0, 

~[—  se^  F2(s)  +  se^TFi(s)]  —  h\D(s)  =  se^F  F4{s)  +  se^  F3(s)]  —  h2D(s) 

—  *1  _  ^1  —  — ^1  _  l±_  — 

e  ci  F2(s)  +  eci  Fi(s)  =  e  c?  F^s)  +  ec2  Fa(s) 

%;[~sF2(s)  +  sF1(s)]  -  h1D{s)  =  0 

se°2  Fi(s)  +  seB2-F,3(s)]  —  h2D(s)  =  0  . 

Since  the  piezoelectric  layers  are  arranged  in  series,  the  total  potential  can  be  written  as  that  for 
capacitors  in  series  via  V (t)  =  V\ (t)  +  V2(t),  where, 


(52) 


Vi (t)  =  f  Ex(x,t)dx 
o 


V2(t)  =  f  Ex(x,t)dx 
h 


(53) 


Substitution  of  equation  26  into  equations  53,  integrating  and  then  Laplace  transforming  the 
result,  allows  determination  of  the  unknown  time-varying  function  D(t),  which  appears  in 


12 


equations  50  as  was  done  in  equation  44.  Prescribing  a  voltage  boundary  condition  V ( t )  allows 
determination  of  the  voltages  via 


Vi(s)  =  —hi [e  ci  F2(s)  +  eci  Fi(s)  -  F2(s)  -  Fi(s)]  +  l1D(s)/e1 

(54) 

V2 (s)  =  -h2[e~^ Fi(s)  +  e^F3(s)  -  e~^~ F4(s)  -  e~^ F3(s)]  +  W(s)/e2  -  liD(s)/e2  . 

Finally,  the  expressions  for  F\ ,  F>,  F>,  F4  and  D(s)  can  be  used  to  obtain  the  solutions  for 
<Ti(x ,  s )  and  <j2(x,  s ).  These  solutions  are  provided  in  appendix  A,  equations  A-l  and  A-2. 


5.  Material  Properties 


From  equations  18  and  19,  it  is  seen  that  eight  material  properties  are  required  for  the  nonlinear 
constitutive  law.  Three  of  these  coefficients  are  the  typical  elastic  moduli,  a;  piezoelectric 
coupling  term,  e;  and  the  dielectric  constant,  e.  Extensive  literature  can  be  found  documenting 
these  material  properties.  For  example,  values  for  alpha-quartz  at  low  temperatures  have  even 
been  measured  (41),  although  the  authors  do  note  that  even  in  the  case  of  quartz,  inconsistent 
values  for  the  elastic  and  piezoelectric  constants  are  often  reported. 

The  remaining  five  coefficients  are  less  readily  available  with  even  more  uncertainty  concerning 
the  accuracy  of  reported  values.  Davison  and  Graham  (42)  present  third-order  elastic  coefficients 
for  several  materials.  However,  they  only  have  third-order  piezoelectric  constants  for  lithium 
niobate,  and  quartz.  These  constants  can  also  be  estimated,  with  appropriate  assumptions,  using 
pressure  derivatives  as  discussed  by  Clayton  (30),  if  such  data  are  available.  While  this  approach 
could  be  used,  for  the  purposes  of  this  report,  it  is  sufficient  to  recognize  that  when  the  response  is 
restricted  to  one  dimension  that  the  distinction  between  different  crystal  classes  is  not  relevant. 
Accordingly,  the  coefficients  of  X-cut  quartz  for  which  these  properties  have  been  experimentally 
determined  have  been  used.  The  properties  are  shown  in  table  3.  The  linear  values  are  also 
shown.  Both  quartz  and  lead  zirconate  titanate  (PZT)  were  used  in  verifying  the  numerical 
implementation  by  comparing  the  predicted  response  with  exact  solutions  to  the  linear  equations. 
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Table  3.  Material  properties  for  quartz  and  PZT. 


Material  Name 

Quartz  (42) 

PZT  (36) 

Present 

D&G 

Value 

Elastic  stiffness  (GPa) 

cE 

riE 

°33 

86.736 

115.41 

Piezoelectric  coupling  (C/m2) 

e 

e33 

0.171 

15.08 

Dielectric  constant  (e°) 

£ 

FV 

6  33 

4.40 

663.2 

Third-order  elastic  (GPa) 

c 

rE 

u333 

-300.0 

Third-order  piezoelectric  (e°  m/V) 

g 

2  ^ijklm 

-1.31 

Electrostrictive  (e°  F/m) 

Q 

/333 

-4.40 

Third-order  dielectric  (F/V) 

V 

fc333 

0(— 3.5  x  10-17) 

Density  (kg/m3) 

p 

2651 

7500 

6.  Results 


The  solution  of  equation  22  was  verified  by  comparing  with  exact  solutions  obtained  by  using  a 
modified  DAC  algorithm  for  numerically  inverting  the  Laplace  transform  (43)  of  the  linear 
piezoelectric  equations  for  a  simple  disk  of  unit  cross-sectional  area  and  unit  thickness. 

Chen  and  Davison  presented  results  for  the  nonlinear  response  of  a  piezoelectric  material  (44). 
They  present  results  where  the  material  properties  remains  elastic  and  nonconducting. 
Consequently,  the  nonlinearity  they  discuss  is  one  of  finite  deformation,  not  material  nonlinearity. 

In  this  report,  we  present  exact  transient  solutions  (in  the  sense  of  a  Laplace  transform  solution) 
for  linear  elastic  piezoelectric  media.  Yang  presents  results  for  a  similar  problem  but  restricted  to 
harmonic  response  (45).  Redwood  also  solves  for  the  transient  response  using  Laplace  transform 
techniques,  though  the  discussion  is  mostly  limited  to  qualitative  results  (46). 

The  poling  direction  was  aligned  with  the  axial  direction  of  the  disk.  In  the  first  verification,  a 
short  circuit  solution  was  obtained  for  a  step  voltage  with  the  boundary  conditions  listed  in  table 
2.  The  voltage  drop  was  1  V  across  the  thickness,  t.  In  the  second  verification,  an  oscillatory 
pressure  boundary  condition  was  applied.  The  frequency,  u,  was  chosen  to  be  close  to  the  axial 
resonance  frequency  of  the  disk. 
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Figures  3  and  4  show  COMSOL  numerical  results  to  a  step  voltage  for  PZT-4  and  quartz, 
respectively.  Figure  3  shows  the  stress  response  at  the  location  x  =  L/2  for  PZT-4.  It  is  seen 
that  the  analytic  solution  obtained  via  the  modified  DAC  techniques  and  the  weak  form  solution 
compare  very  favorably.  Also,  the  response  shown  agrees  qualitatively  with  experimental  results 
obtained  by  Stuetzer  ( 47)  for  PZT-4.  The  effect  of  adding  artificial  damping  is  clearly  evident  as 
the  Rayleigh  damping  completely  eliminates  the  spurious  oscillations.  Similarly,  the  modified 
DAC  algorithm  and  COMSOL  numerical  results  for  the  response  of  a  quartz  disk  subjected  to  a 
unit  step  voltage  are  in  excellent  agreement  (figure  4).  In  the  case  of  quartz,  however,  the 
numerical  oscillations  due  to  the  solver  response  to  step  loadings  is  more  pronounced.  While  the 
inclusion  of  the  numerical  damping  completely  eliminates  these  oscillations,  there  is  a  loss  of 
accuracy  at  the  step  transitions  as  a  larger  value  of  the  damping  parameter  vR  was  required  to 
fully  damp  the  solver  oscillations. 


Undamped 
Rayleigh  Damping 
Modified  DAC 


Figure  3.  Transient  stress  history  (at  x=2.15  mm)  in  a 


4.3-mm-thick  PZT-4  disk  subjected  to  a  Ffeaviside  step 
voltage. 
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Figure  4.  Transient  stress  history  (at  x=2.15  mm)  in  a 

4.3-mm-thick  quartz  disk  subjected  to  a  Heaviside  step 
voltage. 


The  numerical  solution  was  also  verified  by  comparing  the  resonance  response  (figure  5).  As 
sinusoidal  loading  is  reasonably  smooth,  solver  oscillations  are  not  induced.  It  is  seen  that  the 
numerical  and  analytical  solutions  are  again  in  excellent  agreement.  The  effect  of  the  numerical 
damping  is  also  clearly  seen  in  the  reduction  in  the  peak  stresses  of  each  loading  cycle.  Finally, 
the  computed  response  of  a  bi-layered  piezoelectric  media  was  compared  with  the  Laplace 
transform  solution.  The  solution  was  compared  at  two  locations,  1/ 4  (figure  6)  and  3Z/4  (figure 
7),  which  correspond  to  the  midpoints  of  the  PZT  and  quartz  layers,  respectively.  It  is  seen  that 
the  response  is  much  more  complicated  due  to  internal  wave  reflections  due  to  the  impedance 
mismatch  between  the  materials.  As  a  consequence,  the  numerical  solution  exhibits  pronounced 
ringing  effects.  It  is  not  possible  to  reduce  this  ringing  behavior  without  significant  loss  of 
accuracy  in  the  overall  response. 
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Figure  5.  Resonance  response  (at  x=2.15  mm)  in  a  4.3-mm-thick 
quartz  disk  subjected  to  an  harmonic  voltage,  oj  =  666 
kHz. 


Figure  6.  PZT  response  (x=L/4)  in  a  unit  thick  PZT-quartz  disk 
subjected  to  a  step  voltage. 
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Figure  7.  Quartz  response  (x=3L/4)  in  a  unit  thick  PZT-quartz 
disk  subjected  to  a  step  voltage. 


With  the  numerical  solution  technique  verified  by  comparison  with  solutions  obtained  with  the 
modified  DAC  algorithm  solutions,  several  cases  were  studied  corresponding  to  a  step  pressure 
load  applied  at  the  x  =  0  endpoint  location.  The  boundary  conditions  are  given  the  third  column 
of  table  2.  Results  are  shown  in  figure  8.  The  stress  response  is  normalized  in  order  to  highlight 
the  effect  of  pressure  on  the  response.  The  1-GPa  linear  curve  represents  the  exact  solution  of  the 
linear  piezoelectricity  equations  to  a  1-GPa  step  pressure.  The  numerical  solution  of  the 
nonlinear  equations  is  indistinguishable  for  this  load  and  is  therefore  not  shown  in  figure  8.  It  is 
seen  that  the  5-GPa  compression  result  is  also  almost  the  same  as  the  1-GPa  linear  result.  The 
principal  differences  are  that  the  normalized  peak  stress  is  approximately  —2.25  compared  with 
the  peak  value  of  —2.0  for  the  linear  result.  The  wave  speed  is  nearly  the  same  for  these  loads  as 
indicated  by  the  coincidence  of  the  step  changes  in  the  response  as  the  stress  wave  reflects  from 
the  boundaries.  Similarly,  the  5-GPa  tension  and  10-GPa  compression  curves  also  retain  the 
same  basic  form  as  the  linear  response.  However,  the  wave  speed  for  these  two  cases,  while 
different  from  the  linear  response,  are  the  same.  This  indicates  that  the  wave  speed  is  not 
symmetric  with  respect  to  the  loading  level.  Furthermore,  the  wave  speed  is  slower  for  the 
10-GPa  compressive  load  than  the  linear  response.  This  is  in  contrast  with  the  linear  theory, 
which  predicts  that  the  wave  speed  should  increase  with  increasing  compression  load  level.  An 
expression  for  the  effective  wave  speed,  c,  can  be  obtained  by  substituting  equations  18,  19,  and 
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Figure  8.  Transient  stress  histories  (at  x=2.15  mm)  in  a 
4.3-mm-thick  quartz  disk  subjected  to  various 
Heaviside  step  stress  loadings. 


21  into  the  1-D  form  of  equation  of  motion,  equation  8.  This  leads  to 

2  (cE  +  CT  -  gE){  1  +  u,z  )2  +  (aT  -  eE  +  \CT*  -  \QE2  -  gTE) 

c  = - - - - -  .  (55 

Po 

If  terms  arising  in  equation  55  related  to  nonlinear  affects  are  dropped,  then  the  expression  for  c2 
agrees  with  classical  results  (39).  Figure  9  plots  c2  normalized  by  the  linear  elastic  bar  velocity 
versus  the  finite  strain  measure  T.  At  T  =  0,  the  wave  speed  is  the  same  as  classical  elasticity. 
For  slightly  compressive  strains,  T  <  —2.67%,  the  wave  speed  increases.  However,  for  higher 
strain  levels,  the  wave  speed  decreases.  Also,  for  large  finite  strains,  in  either  compression  or 
tension,  the  wave  speed  is  imaginary,  which  could  represent  the  transition  to  standing  waves. 
Whether  this  is  correct  or  simply  represents  a  constraint  on  the  range  of  validity  of  the  theory  is 
currently  being  investigated. 
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Figure  9.  Normalized  wave  speed  squared  c2p/cE  vs.  finite 
strain  Y . 


7.  Conclusions 


A  second-order  theory  for  the  behavior  of  piezoelectric  materials  has  been  presented.  It  has  been 
shown  that  the  presented  theory  reduces  to  the  same  results  as  linear  piezoelectricity  for  small 
strains.  The  effects  of  the  nonlinear  terms  have  been  shown  by  computing  the  response  of  quartz 
to  various  levels  of  pressure,  which  are  realistic  for  shock  loading.  In  these  cases,  the  theory 
predicts  that  for  large  enough  compressive  or  tensile  pressures  that  the  wave  speed  will  decrease. 
This  is  in  contrast  to  linear  piezoelectricity  where  the  wave  speed  always  increases  with 
increasing  compressive  stresses.  Finally,  the  second-order  theory  predicts  that  governing 
equations  will  change  behavior  at  large  absolute  values  of  the  finite  strain.  Whether  this 
represents  a  real  phenomena  or  is  simply  a  restriction  on  the  range  of  applicability  of  the  theory  is 
currently  being  investigated. 
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Appendix  A.  Bi-material  Solution 


The  Laplace  transform  solutions  for  the  stresses  in  each  layer  of  the  bi-material  slab  are  given  by 


<7\{x,s)  =  - 


2Ei£2MiA/l2VoAie  ci  ^C2-/Wi  —  /^2e 2ci  —  c\hiM2A^A\Q^ 


Ai§  +  Ai6  —  A17  +  Ai8 


(A-l) 


a2(x,s)  = 


_S±X±1A  /  _S_  _S_  .  0  A  A  A  \ 

2£i£2MiM2VoAne  c2  [c2h2M\ec^  A^A\2 c\M2e2c2  {J11A2  Ay?,  —  l^A^A^A^ 


A15  +  A16  —  A17  +  A18 


(A-2) 
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where 


Ai  = 

A2  = 
A3  = 
A4  = 
A5  = 

^6  = 
a7  = 
A8  = 

Ag  = 
A\o  = 
An  = 
Al2  = 
A13  — 
Au  = 

Ai5  = 
^16  — 
Air  — 
^4i8  — 


ec2  +  1 
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e^i  —  i 
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ec2  —  1 

>  a 

eci  —  1 

eci  +  1 

>  _S_ 

62c2  _  1 

iX 

eci  —  1 

iX 

e  ci  —  e  2ci 

f  sx 

eci  +  1 


eci  —  eci  ) 
e=2  —  e^2  J 
e=2  —  1 


e  c2  +  e  c2  J 

s  sx  ' 

e  2c2  —  e  c2 


2cl£i£2hlM4  A3A4 

C‘2 M2 A'l i  (4ci£i£2  {h2A^A^  —  hih2A22  A§2  + 

(£i  +  £2)  .MisA.iA.4) 

ciM^As  (2ci£iE2h\A4  —  (£1  +  £2)  MisAs) 


(A- 3) 
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Appendix  B.  Linear  Piezoelectric  Constitutive  Relations 


The  following  forms  of  linear  piezoelectric  constitutive  relationships  are  given  by  the  IEEE 
Standard  on  Peizoelectricity  (48).  There  are  four  standard  forms  of  the  linear  piezoelectric 
constitutive  laws.  Written  in  tensor  notation,  these  forms  are 


Tij  ^-ijkl^kl  &kijDk 

Di  &iklSkl  T  ^ikDk 


(B-l) 


Tij  —  cijki^ki  hkijDk 
Ei  =  —hikiSki  +  /3fkDk 

Sij  =  Sijkl^kl  ~  dkijEk 
Di  dikiTfci  T  c^kEk 


(B-2) 


(B-3) 


Sij  —  sijkl^kl  +  QkijDk 
Ei  =  — gikiTki  +  PjkDk 

Using  Voigt  notation,  the  above  relationships  can  be  written  as 


(B-4) 


Dp  CpqSq  ekpEk 

Di  G{qSq  “t-  ^ik^k 


(B-5) 


Tp  =  CpqSg  ~  hkpDk 
Ei  =  —hiqSq  +  fifkDk 

Sp  =  SpqTq  ~  dkpEk 

Di  =  diqTq  +  eJkEk 

Sp  =  SpqEq  +  QkpDk 
Di  =  —giqTq  +  /3fkDk 

where  the  ranges  on  the  subscripts  are  i,k  =  1,2,3  and  p,  q  =  1, 2,  3, 4,  5, 6. 


(B-6) 


(B-7) 


(B-8) 


The  relationship  between  the  coefficients  is  found  using  equations  B-5-B-8.  In  our  results,  the 
relationship  that  is  required  is  between  the  coefficients  of  the  Stress-Charge  form,  equation  B-5, 
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and  the  Stress- Voltage  form,  equation  B-6. 


Cpq  Cpq  “1“  ekphkp 
kp  h{p6ki 


(B-9) 


In  one  dimension,  this  reduces  to 


cD  =  CE  +  eh 
h=  * 

e 


(B-10) 
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RDRL  WMM  C 
J  LA  SCALA 
RDRL  WMM  D 
ECHIN 
KCHO 

RDRL  WMM  E 
J  ADAMS 
M  COLE 
TIES  SEN 
J  LASALVIA 
P  PATEL 
J  SANDS 
J  SINGH 
RDRL  WMM  F 
L  KECSKES 
H  MAUPIN 
RDRL  WML  G 
I ANDZELM 
A  RAWLETT 
RDRL  WMP 
P  BAKER 
S  SCHOENFELD 


RDRL  WMP  B 
R  BECKER 
S  BILYK 
D  CASEM 
J  CLAYTON 
M  GREENFIELD 
C  HOPPEL 
R  KRAFT 
B  LEAVY 
M  RAFTENBERG 
S  SATAPATHY 
M  SCHEIDLER 
T  WEERASOORIYA 
RDRL  WMP  C 
T  BJERKE 
S  SEGLETES 
RDRL  WMP  D 
R  DONEY 
D  KLEPONIS 
J  RUNYEON 
B  SCOTT 
H  MEYER 
RDRL  WMP  E 
M  BURKINS 
BLOVE 
RDRL  WMP  F 
M  CHOWDHURY 
A  FRYDMAN 
N  GNIAZDOWSKI 
R  GUPTA 
RDRL  WMP  G 
N  ELDREDGE 
D  KOOKER 
S  KUKUCK 
G  R  PEHRSON 
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Intentionally  left  blank. 
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